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[i] Current research has reemphasized the importance of cyclotron resonant wave particle 
interactions for radiation belt electrons. Whistler mode hiss, chorus, and EMIC waves can act in 
combination to cause acceleration and loss of radiation belt electrons at greater rates than 
previously appreciated. These processes can be described by quasi-linear theory, but calculating 
quasi-linear diffusion coefficients is computationally demanding. Recent advances have been 
made in computing bounce averaged quasi-linear pitch angle, energy, and mixed diffusion 
coefficients for hiss and EMIC in the high density plasmasphere; this paper outlines 
generalization of these techniques for chorus waves, prevalent in the low density region outside 
the plasmasphere. These coefficients are associated with a two-dimensional diffusion equation 
whose numerical solution by finite differencing methods requires care, for reasons having to do 
with the relation between the mixed and other diffusion coefficients, as discussed. index terms: 

2716 Magnetospheric Physics: Energetic particles, precipitating; 2720 Magnetospheric Physics: Energetic particles, 
trapped; 7867 Space Plasma Physics: Wave/particle interactions; 2730 Magnetospheric Physics: Magnetosphere— 
inner; KEYWORDS: wave-particle interactions, pitch angle diffusion, diffusion equation, cyclotron resonance 

Citation: Albert, J. M. (2004), Using quasi-linear diffusion to model acceleration and loss from wave-particle interactions. 
Space Weather, 2, S09S03, doi:10.1029/2004SW000069. 


1. Introduction 

[ 2 ] Energetic particles in the inner magnetosphere tend 
to be very stable, because of the three adiabatic invariants 
of their motion [Schulz and Lanzerotti, 1974]. Waves that are 
resonant with one of the three corresponding frequencies 
(gyro, bounce, drift) thus provide an important means of 
affecting the distribution. Quasi-linear theory provides a 
framework for evaluating the effect of waves on test 
particles, assuming that the waves have a broad, continu¬ 
ous distribution in frequency and wavenormal angle. It 
provides diffusion coefficients for pitch angle and energy, 
which apply to the instantaneous position of the particle 
and which may be bounce-averaged. 

[ 3 ] Lyons et al [1972, 1973] applied this approach to 
plasmaspheric hiss acting on radiation belt electrons. They 
showed that the pitch angle diffusion coefficients can 
account for the observed pitch angle distributions. They 
also showed that precipitation lifetimes derived from the 
bounce-averaged pitch angle diffusion rates, combined 
with radial diffusion, reproduces the inner zone/slot/outer 
zone radial structure. These lifetime calculations were 
revisited by Albert [1994] and compared to decay rates 
observed by CRRES [Albert, 2000a]. Recently the calcula¬ 
tions have been redone using more detailed estimates of 
the wave populations originating in hiss, lightning, and 
ground-based VLF transmitters [Abel and Thome, 1998; 
Albert, 1999]. Diffusion by hiss has also been included in 
studies of superthermal electrons [Liemohn et al, 1997], 
electron beams [Khazanov et al, 1999, 2000], and ring 


current electrons [Fok et al, 2001; Khazanov et al, 2003a] 
and protons [Kozyra et al, 1994,1995]. The Salammbo code 
[Beutier and Boscher, 1995; Bourdarie et al, 1996] combines 
pitch angle diffusion from hiss with radial diffusion and 
other sources and losses in a global electron model. 

[ 4 ] Quasi-linear diffusion has also been used to model 
the effects of EMIC (electromagnetic ion cyclotron) waves 
on ring current ions [Jordanova et al, 1996,1997,1998, 2001; 
Khazanov et al, 2002, 2003b] and radiation belt electrons 
[Summers et al, 1998; Summers and Thome, 2003; Albert, 
2003], All of the above calculations treated the waves (hiss 
or EMIC) in the high density approximation appropriate 
inside the plasmasphere. 

[ 5 ] A third important class of waves is known as 

chorus, which propagates in the whistler mode outside 
the plasmasphere. Unlike hiss and EMIC, which mostly 
affect pitch angle, these waves can cause significant diffu¬ 
sion in energy [Summers et al, 1998, 2002; Summers and Ma, 
2000], and may account for the observed rapid regenera¬ 
tion of MeV electrons in the outer zone following their 
depletion by magnetic storms [Meredith et al, 2001, 2002, 
2003]. Energization is expected to be stronger for lower 
ratios of the plasma frequency to gyrofrequency, which is 
large within the plasmasphere but decreases to nearly one 
outside it. Therefore it is important that the most recent 
calculations, which confirm this expectation [Home et al, 
2003], are valid for general values of even though 

results using the high density approximation turn out 
to give surprisingly good agreement (R. Home, personal 
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communication, 2003). Nevertheless, it is desirable to be 
able to evaluate diffusion rates for low density chorus 
properly. An efficient method for doing this is discussed 
in section 3, on the basis of extending techniques devel¬ 
oped for high density hiss. These techniques avoid 
spending computation lime on wavenormal angles and 
harmonic numbers for which the resonant frequency does 
not exist or does not lie between the cutoff frequencies. 
The next logical step is to use the rates in the diffusion 
equation for the particle distribution, but including pitch 
angle, energy, and cross diffusion simultaneously presents 
additional numerical difficulties. Likely reasons for this, 
having to do with the general relation between these 
diffusion coefficients, are discussed in section 4. 

[6] It should be mentioned that quasi-linear theory is not 
the only way to treat cyclotron-resonant wave-particle 
interactions. For narrowband waves, with strongly spatially 
varying parameters, it is more appropriate to consider the 
detailed motion of a particle as it passes through reso¬ 
nance. This approach also has a long history [Roberts and 
Buchsbaum, 1964; Bell, 1984], and a Hamiltonian treatment 
has been given by, e.g., Albert [1993, 2000b, 2002]. The 
relation between the two approaches is discussed by 
Swanson [1989] and by Albert [2001], 

2. Quasi-linear Diffusion 

[ 7 ] The cyclotron resonance condition for a diffusing 
electron is 


w - k\\v y = -nS\/ 7 , 


( 1 ) 


where Vt e = \e\Blmc is the electron's nonrelativistic 
gyrofrequency, y is its relativistic factor, and 0 0 is the 
frequency of the resonant wave. The local pitch angle of 
the particle is a, and the wavenormal angle, between k and 
B, is 0. The parallel wave number k\\ is found from the 
cold plasma index of refraction, r\ = Ar/uj, which is 
determined by uj and 0: 



(RL-PS) sin 2 0 -f 2PS 


± y/(Rr^PSf^V+^D^co^~Q 


/2 PRL. 


( 2 ) 


The standard wave parameters R, L, P, S, and D are 
functions of uj [e.g., Stix, 1962] but do not depend on 0. 
The sign ± is chosen to give the wave mode and 
polarization of interest. 

[8] The local quasi-linear diffusion coefficients are sums 
of integrals over wavenormal angle 0 [Lyons, 1974b]: 

D aa = Dp a = D ap = Dpp = X>; P , (3) 

n n n 


with 


D”a = jdOD n l, D" p = J <?0D£p,DJJp = J dQDf p . (4) 


The integrands are related by 

= D p "°=A 2 D* e a , (5) 

where A n = [-sin acos a/(sin 2 a + nf^/ury)] [Lyons, 1974a]. 
In principle the sums are over an infinite number of n 
values, so a significant amount of computation is required, 
which must be repeated for each value of particle location, 
energy, and pitch angle. Bounce averaging adds another 
layer of integration. 

[ 9 ] Since D"! depends on both 0 and uj, the resonant 
frequency uj(0) must be found from the resonance 
condition in order to carry out the 0 integrations. 
also contains factors to describe the wave magnetic 
field distributions, which are usually modeled as gauss- 
ian in uj and tan 0 with sharp cutoffs: 

f B2 M^(0), 01 < 0 < 02, WLC < W < Wuc 
Cave = < (6) 

[ 0, otherwise 


In general, nothing guarantees in advance that the 
frequency uj(0), once found, will lie between the fre¬ 
quency cutoffs. This can be a drastic source of compu¬ 
tational inefficiency, which can be avoided as described 
below. 

[ 10 ] The resonance condition can be expressed as 


1 

T] 2 


1 ^ 

c 2 (w + wQ e /y) 2 


cos 2 0, 


(7) 


or V(u>, 0) = 0), where from now on V will denote the 

expression on the right-hand side of equation (7) and ^ 
will be used for l/q 2 . For fixed 0, the curve of V versus u; 
has fairly simple behavior, as illustrated in Figure 1. 
With a high density approximation to ^ appropriate 
within the plasmasphere [Lyons, 1974b], the depen¬ 
dence of $ on oj is tractable as well, as discussed by 
Albert [1999] for hiss and by Albert [2003] for EMIC 
waves. Figure 2 indicates the behavior of ^(u) at fixed 0 
for several values of the density parameter uj j; e /fi 2 . 


3. Evaluating the Diffusion Coefficients 

[ 11 ] It is a simple but powerful observation that, at fixed 
0, the curves V and ^ cannot intersect if the smallest V 
value is greater than the largest ^ value, or if the largest 
V value is less than the smallest value. That is, there 
can only be a resonant frequency if 

l^min ^ ^max and V max > (8) 

where the minimum and maximum values are taken 
over the frequency range ujlc to uoltc- The geometric 
behavior of the curves allows one to determine these 
minimum and maximum values. For example, for the 
case n > 0, since V(uj) is always increasing, V^n is V(ujlc) 
and V max is TWc)- Similarly, # min and ^ max can be 
found explicitly, using ^ appropriate to high density 
hiss and EMIC waves. It then turns out that conditions 
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Figure 1. The behavior of the function V(w), indicated by plotting VIV 0 against l o/Q e , where Vo is 
(u^/c^cos 2 © (see equation 7). The solid lines are for positive values of n and the dashed lines are 
for negative values of n. Curves with negative n have a singularity at uj/Q e = \n\ly. For large w, the 
curves all approach the n = 0 value V = V 0 (dotted line). Here, for illustration, the relativistic factor 
y is 5, corresponding to about a 2 MeV electron. 


like those in (8) can be inverted algebraically to give 
inequalities for cos 2 0 of the form 

A cos 4 0 + B cos 2 0 -f C > 0, (9) 

where A, B, and C are complicated but fully specified 
functions of uj lc and w U c- The corresponding 0 values 
are easily found by considering the zeroes of the 
quadratic and the sign of A. Only these values will be 
compatible with the conditions in equation (8); all other 
0 ranges can be omitted from the integrals for the 
diffusion coefficients. 

[ 12 ] Since V decreases as \n\ increases, and T' is inde¬ 
pendent of n, there will be a value of \n\ large enough 
that V max is less than \Fmi n for all 0. This and all larger 
values of \n\ will not contribute any resonances to D at 
all. This gives a systematic criterion for cutting off the 
infinite sum over n in equation (3). 

[ 13 ] The success of this approach hinges on being able to 
characterize the behavior of ^(w), so that ^ max and ^ m in 
can be specified for any general value of 0. As men¬ 
tioned above, this has previously been carried out for 
plasmaspheric hiss and EMIC waves. For whistler 
mode chorus with general WpJQe, ^(u) has the same 
qualitative shape as in the high density case; this has 
been verified through extensive numerical evaluation. 
Therefore \F min is the smaller of \F(w LC ) and ^(wucX but 
^max = ^iupeak), which can be found exactly in the high 
density limit, can only be estimated. One such estimate 
is the high density value, because \F(io) < \F HD (ou) so that 



0.0 0.2 0.4 0.6 0.8 1.0 

w/Q e 


Figure 2. The behavior of the function ^(oj) = l/rj 2 for 
fixed 0, indicated by plotting ( 10 2 /Q 2 )# against uj IQ, e at 
0 = 20° for several values of u )p e l&e (see equation 2). The 
top curve shows the high density approximation 
appropriate to whistler mode waves within the plasma- 
sphere. The cyclotron resonance condition can be 
expressed as V(oj, 0) = \F(u>, 0). 
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^max < peak )• More refined upper bounds can be 

developed, and will be reported in detail elsewhere. 
With these estimates of and \F max , the 0 ranges in 
(4) may be reduced, and the diffusion coefficients 
calculated. This technique does not change the results 
at all, but allows them to be found much more effi¬ 
ciently, since computation time is not spent trying to 
obtain resonant frequencies which do not exist or 
which do not lie between the cutoff values, nor spent 
considering unnecessarily large values of ±n. 

4. Diffusion Coefficients and 
Diffusion Equations 

[ 14 ] If n ^ 0 and wy < 0,# the definition of A„ and 
equation (5) show that D”° a > |D”®| > D"®. Equation (5) 
also implies 


D aa D pp = ( D Sp) 2 - ( 10 ) 


[ 16 ] The bounce-averaged version of the diffusion equa¬ 
tion is [e.g., Kozyra et al, 1994] 


Of 


£ = JL A G fn 1 

dt Gpdom \ aoao p dao 

, IA G fn 1 it, D d l 

Gdp \ pa °pdoio pp dp 


(15) 


where G = p 2 T( ao)sin a 0 cos Oq [Schulz, 1991] and T(ao) « 
1.30-0.56 sin olq is the normalized bounce period [e.g., 
Lyons et al., 1972]. The bounce averaged diffusion coeffi¬ 
cients are 


60*00 = J d «o(^) 2 * = J D aa W(\)d\, (16) 

etc., where W(\) = cos a cos 7 X/T(ao)cos 2 olq is the bounce 
average weighting factor [Lyons et al, 1972]. The bounce- 
average equivalent of (11) is 


On the other hand, ignoring variations in the coeffi¬ 
cients, the condition 


AxoaoDpp > 



(17) 


Dolol Dpp > 


( 11 ) 


must be satisfied in order for the 2D diffusion equation 
[Lyons, 1974b] 


d l 

dt 


sinaf D, 


1V +D df 

p da ap dp 


p sin a da 

D 1 y +D 

p 2 dpP \ p(X pda pp dp) 


+ 


( 12 ) 


to be numerically well-behaved [Richtmyer and Morton, 
1967, Gourlay and McKee, 1977]. 

[ 15 ] In general, the integral version of the Cauchy- 
Schwarz inequality is 


JfW* [|^(e)de| > [//(e),(6 } J 2 


(13) 


for any functions/and g, with equality holding only if/is a 
constant times g [Gradshteyn and Ryzhik, 1980]. Taking/(0) = 
(D”„) 1/2 and g(0) = (DSL)*® A n gives 


D n D n > 

on pp ^ 

K) 2 . 

n/0, 

D n D n = 

aa pp 

K) 2 . 

n = 0, 


since A n depends on 0 through 00 except when n - 0. 
Summing over n, (11) is satisfied as long as resonances 
with n^0 are present. Otherwise, the diffusion is not 
truly two dimensional, and equation (12) can be 
replaced by a ID diffusion equation in p\\ [Lyons et al, 
1972]. 


For a single n, a version of the Cauchy-Schwarz inequality 
can be applied to integrals over X, with f{\) = (D act W) 1/2 
and g(\) = ( D pp W ) 112 . Since W depends on X, (17) is 
confirmed even for n = 0. However, n = 0 resonances 
typically occur near the mirror point, where v\\!c = l/rj|j is 
small, not over a wide range of latitude. Thus if only n- 0 
resonances are present, the bounce-averaged diffusion 
coefficients are almost the same as their local, mirror point 
values and (17) is just barely satisfied, so again the 
diffusion is nearly one-dimensional (in the second adia¬ 
batic invariant J). 

[ 17 ] Transforming variables from ( 0 ^, p) to some other 
set (u, v), such as the adiabatic invariants (\l, ]), does not 
change things since 


D ttM D w — D? — 


d(u,v) 


d(oo ,p) 


(&c toaoDpp 


(18) 


so that for nonsingular transformations the 2D condition is 
the same in either set of variables. As a consequence pure 
pitch angle diffusion, which violates (17), is problematic 
when treated as two-dimensional in (jx, J). The Salammbo 
code, which took this approach, only ran stably when 
strong radial diffusion was also included (S. Bourdarie, 
personal communication, 2002) though it has recently been 
rewritten in (a 0/ p) [Boscher, 2003]. Kozyra et al [1995] were 
led to rewrite equation (15) (for nonrelativistic ions) in 
conservative form and use a finite volume scheme, but the 
results are inconclusive because of coding errors in 
evaluating the diffusion coefficients [Liemohn et al, 1997]. 
Other models have typically treated diffusion as ID in oq 
\jordanova et al, 1997; Khazanov et al, 2003b], though some 
simulations included energy diffusion without the cross 
term [Jordanova et al, 1998, 2001], One current effort recasts 
the 2D diffusion process in (\i, J) as a nonlinear advection 
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problem and relies on a sophisticated flux-limiting algo¬ 
rithm (M.-C. Fok, personal communication, 2004), but no 
results have yet been published. Thus the reliable 
numerical solution of equation (15) or its equivalent 
remains an open problem. 

[is] Treating the diffusion coefficients as constant or 
slowly varying, stability analysis of standard finite 
differencing methods [Richtmyer and Morton , 1967] shows 
that the fully explicit scheme should be stable for timestep 
At small enough that 


D u 


r + ■ 


(Au) 2 (At;) 2 


*4 


(19) 


and that the fully implicit and Crank-Nicholson schemes 
should be unconditionally stable for any choice of 
variables (w, v). Nevertheless one expects, and experience 
suggests, that numerical behavior is better with variables 
that make the "cross term" D uv and its gradients small 
compared to the larger of the "diagonal terms" D uu and 
D vv . From this point of view, (olq, p) is probably preferable 
to (p, J) [Kamey, 1986; Reveille et al, 2001]. 


5. Summary 

[ 19 ] As a description of wave-particle interactions, quasi- 
linear diffusion is of renewed importance in understand¬ 
ing radiation belt dynamics. The techniques described 
here make it feasible to compute the diffusion coefficients 
for low density chorus waves, as well as whistler mode hiss 
and broadband EMIC waves, by identifying wavenormal 
angle ranges of waves within a prescribed frequency band 
that are or are not resonant with specified particles. These 
coefficients are for diffusion in equatorial pitch angle c% 
momentum p, and mixed diffusion (D^p) and are to be 
used to advance the time-dependent diffusion equation, 
but this is a surprisingly challenging problem in its own 
right. It was argued above that a source of numerical 
problems is that, because of the relation (5) between the 
different diffusion coefficients, the underlying character of 
the diffusion can actually be one-dimensional or nearly so; 
also, working in the variables (|x, ]) rather than (olq, p) is 
likely to make the numerical stability problems worse. 
This is unfortunate since (jx, J) is certainly more convenient 
for including radial diffusion at constant \x and J. 

[ 20 ] Acknowledgments. This work was supported by the Space 
Vehicles Directorate of the Air Force Research Laboratory and by the 
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